library(ggplot2)
library(reshape2)
library(dplyr)
library(tidyr)
library(GGally)
library(grid)
"%&%" = function(a,b) paste(a,b,sep="")
se <- function(x) sqrt(var(x,na.rm=TRUE)/length(is.na(x)==FALSE))
source('/Volumes/im-lab/nas40t2/hwheeler/PrediXcan_CV/GTEx_2014-06013_release/transfers/PrediXmod/Paper_plots/multiplot.R')
my.dir <- '/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/'
fig.dir <- '~/GitHub/GenArch/GenArchPaper/Figures/'
rna.dir <- my.dir %&% "gtex-rnaseq/"
annot.dir <- my.dir %&% "gtex-annot/"
out.dir <- rna.dir %&% "ind-tissues-RPKM/"
h2.dir <- my.dir %&% "gtex-h2-estimates/"
tislist <- scan(my.dir %&% '40.tissue.list',sep="\n",what="character")
table1 <- matrix(0,nrow=length(tislist)+1,ncol=6)
ct <- read.table(my.dir %&% 'cross-tissue.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T)
n <- ct$N[1]
numexpgenes <- dim(ct)[1]
numexpgenes <- length(ct$local.p[is.na(ct$local.p)==FALSE])
meanh2 <- sprintf("%.3f",mean(ct$local.h2,na.rm=TRUE))
semean <- sprintf("%.4f",se(ct$local.h2))
meanandse <- meanh2 %&% " (" %&% semean %&% ")"
pest <- ct %>% mutate(pchi=pchisq((local.h2/local.se)^2, df=1, lower.tail=FALSE)) %>% mutate(local.q=p.adjust(pchi,method="BH")) %>% arrange(local.h2) %>% mutate(Qlt05=local.q<0.1)
propsig <- sprintf("%.1f",table(pest$Qlt05)[2]/sum(table(pest$Qlt05,useNA="i"))*100)
numsig <- table(pest$Qlt05)[2]
table1[1,] <- c("Cross-tissue",n,meanh2,propsig,numsig,numexpgenes)
hist(pest$local.p,main="CT")

hist(pest$pchi,main="CT")

for(i in 1:length(tislist)){
tis <- tislist[i]
data <- read.table(my.dir %&% 'GTEx.TS.' %&% tis %&% '.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T,sep="\t")
n <- data$N[1]
numexpgenes <- dim(data)[1]
numexpgenes <- length(data$local.p[is.na(data$local.p)==FALSE]) ##num expressed genes mean(RPKM)>0.1
meanh2 <- sprintf("%.3f",mean(data$local.h2,na.rm=TRUE))
semean <- sprintf("%.4f",se(data$local.h2))
pest <- data %>% mutate(pchi=pchisq((local.h2/local.se)^2, df=1, lower.tail=FALSE)) %>% mutate(local.q=p.adjust(pchi,method="BH")) %>% arrange(local.h2) %>% mutate(Qlt05=local.q<0.1)
propsig <- sprintf("%.2f",table(pest$Qlt05)[2]/sum(table(pest$Qlt05,useNA="i"))*100)
numsig <- table(pest$Qlt05)[2]
meanandse <- meanh2 %&% " (" %&% semean %&% ")"
tableinfo <- c(tis,n,meanh2,propsig,numsig,numexpgenes)
table1[i+1,] <- tableinfo
hist(pest$local.p,main=tis)
hist(pest$pchi,main=tis)
}
















































































colnames(table1)=c("tissue","n","mean h2","% FDR<0.1","num FDR<0.1","num expressed")
#table1
library(xtable)
tab <- xtable(table1)
print(tab, type="latex",include.rownames=FALSE)
## % latex table generated in R 3.2.0 by xtable 1.8-0 package
## % Wed Aug 10 21:44:04 2016
## \begin{table}[ht]
## \centering
## \begin{tabular}{llllll}
## \hline
## tissue & n & mean h2 & \% FDR$<$0.1 & num FDR$<$0.1 & num expressed \\
## \hline
## Cross-tissue & 450 & 0.062 & 17.7 & 2995 & 14861 \\
## Adipose - Subcutaneous & 298 & 0.017 & 8.31 & 1408 & 14132 \\
## Adrenal Gland & 126 & 0.037 & 7.84 & 1329 & 14268 \\
## Artery - Aorta & 198 & 0.022 & 8.55 & 1449 & 14634 \\
## Artery - Coronary & 119 & 0.048 & 7.14 & 1211 & 14663 \\
## Artery - Tibial & 285 & 0.022 & 7.44 & 1262 & 14534 \\
## Brain - Anterior cingulate cortex (BA24) & 72 & 0.037 & 9.60 & 1627 & 14904 \\
## Brain - Caudate (basal ganglia) & 100 & 0.042 & 8.19 & 1388 & 14562 \\
## Brain - Cerebellar Hemisphere & 89 & 0.046 & 9.60 & 1627 & 14851 \\
## Brain - Cerebellum & 103 & 0.041 & 9.55 & 1618 & 14675 \\
## Brain - Cortex & 96 & 0.053 & 8.73 & 1480 & 14405 \\
## Brain - Frontal Cortex (BA9) & 92 & 0.037 & 9.26 & 1569 & 14481 \\
## Brain - Hippocampus & 81 & 0.033 & 9.55 & 1618 & 14504 \\
## Brain - Hypothalamus & 81 & 0.020 & 10.85 & 1839 & 14563 \\
## Brain - Nucleus accumbens (basal ganglia) & 93 & 0.040 & 9.35 & 1585 & 14546 \\
## Brain - Putamen (basal ganglia) & 82 & 0.033 & 9.33 & 1581 & 14319 \\
## Breast - Mammary Tissue & 183 & 0.020 & 8.06 & 1367 & 14080 \\
## Cells - EBV-transformed lymphocytes & 115 & 0.045 & 7.80 & 1323 & 14072 \\
## Cells - Transformed fibroblasts & 272 & 0.019 & 7.87 & 1334 & 14193 \\
## Colon - Sigmoid & 124 & 0.024 & 9.31 & 1578 & 14749 \\
## Colon - Transverse & 170 & 0.022 & 9.28 & 1573 & 14217 \\
## Esophagus - Gastroesophageal Junction & 127 & 0.029 & 8.01 & 1358 & 14581 \\
## Esophagus - Mucosa & 242 & 0.026 & 7.20 & 1220 & 14944 \\
## Esophagus - Muscularis & 219 & 0.024 & 7.93 & 1344 & 14805 \\
## Heart - Atrial Appendage & 159 & 0.031 & 8.48 & 1437 & 14401 \\
## Heart - Left Ventricle & 190 & 0.025 & 7.82 & 1325 & 14703 \\
## Liver & 98 & 0.036 & 8.96 & 1519 & 14476 \\
## Lung & 279 & 0.018 & 7.72 & 1308 & 14452 \\
## Muscle - Skeletal & 361 & 0.020 & 7.92 & 1342 & 14382 \\
## Nerve - Tibial & 256 & 0.026 & 7.29 & 1235 & 14347 \\
## Ovary & 85 & 0.043 & 8.98 & 1523 & 13922 \\
## Pancreas & 150 & 0.036 & 8.24 & 1396 & 14678 \\
## Pituitary & 87 & 0.044 & 8.74 & 1481 & 14752 \\
## Skin - Not Sun Exposed (Suprapubic) & 196 & 0.041 & 5.53 & 938 & 15062 \\
## Skin - Sun Exposed (Lower leg) & 303 & 0.027 & 6.02 & 1020 & 14808 \\
## Small Intestine - Terminal Ileum & 77 & 0.046 & 10.22 & 1733 & 13820 \\
## Spleen & 89 & 0.062 & 8.54 & 1448 & 14230 \\
## Stomach & 171 & 0.020 & 9.04 & 1532 & 14259 \\
## Testis & 157 & 0.042 & 9.21 & 1561 & 14779 \\
## Thyroid & 279 & 0.024 & 7.63 & 1293 & 14478 \\
## Whole Blood & 339 & 0.025 & 6.78 & 1150 & 14609 \\
## \hline
## \end{tabular}
## \end{table}